Equation Solving

One-dimensional solvers

Solution of one equation in one variable is accomplished by children of the class o2scl::root.

For one-dimensional solving, if the root is bracketed, use o2scl::root_bkt_cern or o2scl::root_brent_gsl. The o2scl::root_bkt_cern class is typically faster (for the same accuracy) than o2scl::root_brent_gsl. If a relatively fast derivative is available, use o2scl::root_stef. If neither a bracket nor a derivative is available, you can use o2scl::root_cern.

The o2scl::root base class provides the structure for three different solving methods:

  • o2scl::root::solve() which solves a function given an initial guess x
  • o2scl::root::solve_bkt() which solves a function given a solution bracketed between x1 and x2. The values of the function at x1 and x2 should have different signs.
  • o2scl::root::solve_de() which solves a function given an initial guess x and the function's derivative.

There is an example using the one-dimensional solver in the Function object example .

The o2scl::root base class also contains the relative tolerance (o2scl::root::tol_rel), absolute tolerance (o2scl::root::tol_abs), the number of iterations (o2scl::root::ntrial), the verbosity parameter (o2scl::root::verbose), and the number of iterations in the last solve (o2scl::root::last_ntrial).

If not all of these three functions are overloaded, then the source code in the o2scl::root base class is designed to try to automatically provide the solution using the remaining functions. Most of the one-dimensional solving routines, in their original form, are written in the second or third form above. For example, o2scl::root_brent_gsl is originally a bracketing routine of the form o2scl::root::solve_bkt(), but calls to either o2scl::root::solve() or o2scl::root::solve_de() will attempt to automatically bracket the function given the initial guess that is provided. Of course, it is frequently most efficient to use the solver in the way it was intended.

Multi-dimensional solvers

Solution of more than one equation is accomplished by descendants of the class o2scl::mroot . The higher-level interface is provided by the function o2scl::mroot::msolve() .

For multi-dimensional solving, you can use either o2scl::mroot_cern or o2scl::mroot_hybrids. While o2scl::mroot_cern cannot utilize user-supplied derivatives, o2scl::mroot_hybrids can use user-supplied derivative information (as in the GSL hybridsj method) using the function o2scl::mroot_hybrids::msolve_de() .

A specialization of o2scl::mroot_hybrids for Armadillo is given in o2scl::mroot_hybrids_arma_qr_econ where the QR decomposition used in the solver is performed by the Armadillo library. A similar specialization for Eigen is in o2scl::mroot_hybrids_eigen . These specializations will be faster than when the number of variables is sufficiently large.

Multi-dimensional solver example

This demonstrates several ways of using the multi-dimensional solvers to solve the equations

\begin{eqnarray*} \sin \left( x_0 - \frac{1}{4} \right) &=& 0 \nonumber \\ \sin \left( x_1 - \frac{1}{5} \right) &=& 0 \end{eqnarray*}

/* Example: ex_mroot.cpp
-------------------------------------------------------------------
Several ways to use an O2scl solver to solve a simple function
*/
#include <cmath>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <o2scl/test_mgr.h>
#include <o2scl/mm_funct.h>
#include <o2scl/mroot_hybrids.h>
#include <o2scl/mroot_cern.h>
using namespace std;
using namespace o2scl;
int gfn(size_t nv, const ubvector &x, ubvector &y) {
y[0]=sin(x[1]-0.2);
y[1]=sin(x[0]-0.25);
return 0;
}
class cl {
public:
// Store the number of function and derivative evaluations
int nf, nd;
int mfn(size_t nv, const ubvector &x, ubvector &y) {
y[0]=sin(x[1]-0.2);
y[1]=sin(x[0]-0.25);
nf++;
return 0;
}
int operator()(size_t nv, const ubvector &x, ubvector &y) {
y[0]=sin(x[1]-0.2);
y[1]=sin(x[0]-0.25);
nf++;
return 0;
}
int mfnd(size_t nx, ubvector &x, size_t ny,
ubvector &y, ubmatrix &j) {
j(0,0)=0.0;
j(0,1)=cos(x[1]-0.2);
j(1,0)=cos(x[0]-0.25);
j(1,1)=0.0;
nd++;
return 0;
}
template<class vec_t>
int mfn_tlate(size_t nv, const vec_t &x, vec_t &y) {
y[0]=sin(x[1]-0.2);
y[1]=sin(x[0]-0.25);
nf++;
return 0;
}
template<class vec_t, class mat_t>
int mfnd_tlate(size_t nx, vec_t &x, size_t ny, vec_t &y, mat_t &j) {
j(0,0)=0.0;
j(0,1)=cos(x[1]-0.2);
j(1,0)=cos(x[0]-0.25);
j(1,1)=0.0;
nd++;
return 0;
}
};
int main(void) {
cout.setf(ios::scientific);
cl acl;
ubvector x(2);
/*
Using a member function with ublas vectors
*/
mm_funct f1=std::bind
(std::mem_fn<int(size_t,const ubvector &,ubvector &)>(&cl::mfn),
&acl,std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
x[0]=0.5;
x[1]=0.5;
acl.nf=0;
int ret1=cr1.msolve(2,x,f1);
cout << "GSL solver (numerical Jacobian): " << endl;
cout << "Return value: " << ret1 << endl;
cout << "Number of iterations: " << cr1.last_ntrial << endl;
cout << "Number of function evaluations: " << acl.nf << endl;
cout << endl;
t.test_rel(x[0],0.25,1.0e-6,"1a");
t.test_rel(x[1],0.2,1.0e-6,"1b");
/*
Using the CERNLIB solver
*/
x[0]=0.5;
x[1]=0.5;
acl.nf=0;
int ret2=cr2.msolve(2,x,f1);
cout << "CERNLIB solver (numerical Jacobian): " << endl;
cout << "Return value: " << ret2 << endl;
cout << "INFO parameter: " << cr2.get_info() << endl;
cout << "Number of function evaluations: " << acl.nf << endl;
cout << endl;
t.test_rel(x[0],0.25,1.0e-6,"2a");
t.test_rel(x[1],0.2,1.0e-6,"2b");
/*
Using a member function with \ref ovector objects, but
using the GSL-like interface with set() and iterate().
*/
x[0]=0.5;
x[1]=0.5;
cr3.allocate(2);
cr3.set(2,x,f1);
bool done=false;
do {
double r3=cr3.iterate();
double resid=fabs(cr3.f[0])+fabs(cr3.f[1]);
if (resid<cr3.tol_rel || r3>0) done=true;
} while (done==false);
t.test_rel(cr3.x[0],0.25,1.0e-6,"3a");
t.test_rel(cr3.x[1],0.2,1.0e-6,"3b");
/*
Now instead of using the automatic Jacobian, using
a user-specified Jacobian.
*/
jac_funct j4=std::bind
(std::mem_fn<int(size_t,ubvector &,size_t,ubvector &,ubmatrix &)>
(&cl::mfnd),&acl,std::placeholders::_1,std::placeholders::_2,
std::placeholders::_3,std::placeholders::_4,std::placeholders::_5);
x[0]=0.5;
x[1]=0.5;
acl.nf=0;
acl.nd=0;
int ret4=cr1.msolve_de(2,x,f1,j4);
cout << "GSL solver (analytic Jacobian): " << endl;
cout << "Return value: " << ret4 << endl;
cout << "Number of iterations: " << cr1.last_ntrial << endl;
cout << "Number of function evaluations: " << acl.nf << endl;
cout << "Number of Jacobian evaluations: " << acl.nd << endl;
cout << endl;
t.test_rel(x[0],0.25,1.0e-6,"4a");
t.test_rel(x[1],0.2,1.0e-6,"4b");
/*
Using a user-specified Jacobian and the GSL-like interface
*/
x[0]=0.5;
x[1]=0.5;
cr5.allocate(2);
cr5.set_de(2,x,f1,j4);
done=false;
do {
double r3=cr5.iterate();
double resid=fabs(cr5.f[0])+fabs(cr5.f[1]);
if (resid<cr5.tol_rel || r3>0) done=true;
} while (done==false);
t.test_rel(cr5.x[0],0.25,1.0e-6,"5a");
t.test_rel(cr5.x[1],0.2,1.0e-6,"5b");
/*
Using a class with an operator(). Note that there can be only one
operator() function in each class.
*/
x[0]=0.5;
x[1]=0.5;
cr9.msolve(2,x,acl);
t.test_rel(x[0],0.25,1.0e-6,"9a");
t.test_rel(x[1],0.2,1.0e-6,"9b");
/*
Using a function pointer to a global function.
*/
typedef int (*gfnt)(size_t, const ubvector &, ubvector &);
gfnt f10=&gfn;
x[0]=0.5;
x[1]=0.5;
cr10.msolve(2,x,f10);
t.test_rel(x[0],0.25,1.0e-6,"10a");
t.test_rel(x[1],0.2,1.0e-6,"10b");
/*
Using std::vector<double>
*/
std::vector<double> svx(2);
svx[0]=0.5;
svx[1]=0.5;
typedef std::function<int(size_t,const std::vector<double> &,
std::vector<double> &) > mm_funct_sv;
typedef std::function<
int(size_t,std::vector<double> &,size_t,std::vector<double> &,
ubmatrix &) > jac_funct_sv;
mm_funct_sv f11=std::bind
(std::mem_fn<int(size_t,const std::vector<double> &,
std::vector<double> &)>
(&cl::mfn_tlate<std::vector<double> >),
&acl,std::placeholders::_1,std::placeholders::_2,std::placeholders::_3);
ubmatrix,jac_funct_sv> cr11;
cr11.msolve(2,svx,f11);
t.test_rel(x[0],0.25,1.0e-6,"11a");
t.test_rel(x[1],0.2,1.0e-6,"11b");
/*
Using ublas::unbounded_array<double>
*/
typedef boost::numeric::ublas::unbounded_array<double> uarr;
typedef std::function<int(size_t,const uarr &,uarr &)> mm_funct_ua;
typedef std::function<int(size_t,uarr &,size_t,uarr &,
ubmatrix &) > jac_funct_ua;
uarr uad(2);
uad[0]=0.5;
uad[1]=0.5;
mm_funct_ua f12=std::bind(std::mem_fn<int(size_t,const uarr &,uarr &)>
(&cl::mfn_tlate<uarr>),
&acl,std::placeholders::_1,
std::placeholders::_2,std::placeholders::_3);
cr12.msolve(2,uad,f12);
t.test_rel(x[0],0.25,1.0e-6,"12a");
t.test_rel(x[1],0.2,1.0e-6,"12b");
/*
Using ublas::bounded_array<double>
*/
typedef boost::numeric::ublas::bounded_array<double,2> barr;
typedef std::function<int(size_t,const barr &,barr &)> mm_funct_ba;
typedef std::function<int(size_t,barr &,size_t,barr &,
ubmatrix &) > jac_funct_ba;
barr bad(2);
bad[0]=0.5;
bad[1]=0.5;
mm_funct_ba f13=std::bind(std::mem_fn<int(size_t,const barr &,barr &)>
(&cl::mfn_tlate<barr>),
&acl,std::placeholders::_1,
std::placeholders::_2,std::placeholders::_3);
cr13.msolve(2,bad,f13);
t.test_rel(x[0],0.25,1.0e-6,"13a");
t.test_rel(x[1],0.2,1.0e-6,"13b");
#ifdef O2SCL_EIGEN
typedef std::function<int(size_t, const Eigen::VectorXd &,
Eigen::VectorXd &)> mm_funct_eigen;
typedef std::function<int(size_t,Eigen::VectorXd &,size_t,Eigen::VectorXd &,
ubmatrix &) > jac_funct_eigen;
// Using Eigen with an analytic Jacobian
mm_funct_eigen f14=std::bind
(std::mem_fn<int(size_t,const Eigen::VectorXd &,Eigen::VectorXd &)>
(&cl::mfn_tlate<Eigen::VectorXd>),
&acl,std::placeholders::_1,
std::placeholders::_2,std::placeholders::_3);
jac_funct_eigen fd14=std::bind
(std::mem_fn<int(size_t,Eigen::VectorXd &,size_t,
Eigen::VectorXd &,Eigen::MatrixXd &)>
(&cl::mfnd_tlate<Eigen::VectorXd>),&acl,
std::placeholders::_1,std::placeholders::_2,
std::placeholders::_3,std::placeholders::_4,std::placeholders::_5);
mroot_hybrids<mm_funct_eigen,Eigen::VectorXd,
Eigen::MatrixXd,jac_funct_eigen> cr14;
Eigen::VectorXd xE(2);
xE[0]=0.5;
xE[1]=0.5;
cr14.msolve_de(2,xE,f14,fd14);
t.test_rel(xE[0],0.25,1.0e-6,"eigen a");
t.test_rel(xE[1],0.2,1.0e-6,"eigen b");
#endif
t.report();
return 0;
}
// End of example
boost::numeric::ublas::matrix< double >
o2scl::mroot_cern
Multi-dimensional mroot-finding routine (CERNLIB)
Definition: mroot_cern.h:95
boost::numeric::ublas::vector< double >
o2scl::mroot_hybrids::x
vec_t x
The present solution.
Definition: mroot_hybrids.h:697
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::test_mgr::report
bool report() const
Provide a report of all tests so far.
o2scl::mroot_hybrids::iterate
int iterate()
Perform an iteration.
Definition: mroot_hybrids.h:704
o2scl::mroot_cern::get_info
int get_info()
Get the value of INFO from the last call to msolve()
Definition: mroot_cern.h:166
o2scl::mroot< mm_funct, boost::numeric::ublas::vector< double >, jac_funct >::last_ntrial
int last_ntrial
The number of iterations for in the most recent minimization.
Definition: mroot.h:94
o2scl::mroot_hybrids::f
vec_t f
The value of the function at the present iteration.
Definition: mroot_hybrids.h:694
o2scl::test_mgr::test_rel
bool test_rel(data_t result, data_t expected, data_t rel_error, std::string description)
Test for .
Definition: test_mgr.h:168
o2scl::mroot_hybrids::msolve
virtual int msolve(size_t nn, vec_t &xx, func_t &ufunc)
Solve ufunc using xx as an initial guess, returning xx.
Definition: mroot_hybrids.h:996
o2scl::mroot_hybrids::allocate
void allocate(size_t n)
Allocate the memory.
Definition: mroot_hybrids.h:922
o2scl::mroot_hybrids
Multidimensional root-finding algorithm using Powell's Hybrid method (GSL)
Definition: mroot_hybrids.h:521
o2scl::test_mgr
A class to manage testing and record success and failure.
Definition: test_mgr.h:50
o2scl::mroot_hybrids::set_de
int set_de(size_t nn, vec_t &ax, func_t &ufunc, jfunc_t &dfunc)
Set the function, the Jacobian, the parameters, and the initial guess.
Definition: mroot_hybrids.h:1078
o2scl::mroot_cern::msolve
virtual int msolve(size_t nvar, vec_t &x, func_t &func)
Solve func using x as an initial guess, returning x.
Definition: mroot_cern.h:233
o2scl::mroot_hybrids::msolve_de
virtual int msolve_de(size_t nn, vec_t &xx, func_t &ufunc, jfunc_t &dfunc)
Solve func with derivatives dfunc using x as an initial guess, returning x.
Definition: mroot_hybrids.h:983
o2scl::mm_funct
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef in src/base/mm_funct.h.
Definition: mm_funct.h:43
o2scl::test_mgr::set_output_level
void set_output_level(int l)
Set the output level.
Definition: test_mgr.h:119
o2scl::jac_funct
std::function< int(size_t, boost::numeric::ublas::vector< double > &, size_t, boost::numeric::ublas::vector< double > &, boost::numeric::ublas::matrix< double > &) > jac_funct
Jacobian function (not necessarily square) in src/root/jacobian.h.
Definition: jacobian.h:44
o2scl::mroot_hybrids::set
int set(size_t nn, vec_t &ax, func_t &ufunc)
Set the function, the parameters, and the initial guess.
Definition: mroot_hybrids.h:1011

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).