One Dimensional Root-Finding

Feature

Currently I am developing the routines for finding roots of arbitrary one-dimensional functions. The idea is to offer simple functions that allow it, without the need to implement a robust search system such as workspaces and a little bigger structures like thoose.

Discussion

Since I do not have much knowledge about the different existing methods for the search of roots, would be grateful for any type of contribution that could be made.

In the first comment I show the methods that I plan to implement in this new feature as well as the prototype of each of them with the idea that each one of them is up for discussion.

Contributing

Those who wish to contribute to this module are asked to discuss the implementation of existing methods as well as the suggestion and development of new ones.

cml/roots.h

[. . .]

double cml_root_brent(cml_function_t *func, double x1, double x2, double *tol);

int cml_root_newton_bisection(cml_function_fdf_t *func, double x_min, double x_max,
                              double tol,int n_max,double * res);

int cml_root_newton(cml_function_fdf_t *func, double x0,
                    double epsrel, double epsabs,
                    int n_max, double * res);

int cml_root_bisection(cml_function_t *funcunc, double xmin, double xmax,
                       double epsrel, double epsabs,
                       size_t n_max, double *res);

int cml_multiroot_newton(cml_function_vec_t *funcunc, const double *x0,
                         double x_eps, double fx_eps, int n_max,
                         bool verbose, double *res);

int cml_root_fsolve(cml_function_vec_t *func, double *x, double *fx,
                    double xtol, int maxfev, int *nfev,
                    double *scale, int error_msg);

int cml_root_fsolve_lsq(cml_function_vec_t *func, double *x, int m,
                        double *fx,  double xtol, double ftol,
                        double gtol, int maxfev, int *nfev,
                        double *scale, int error_msg);

[. . .]

Currently implemented,

/**
 * Find the root of a function using a bisection method
 *
 * @param xmin lower bound
 * @param xmax upper bound
 * @param epsrel if the relative improvement over the root is less than this value,
 * then break;
 * @param epsabs if the absolute improvement over the root is less than this value,
 * then break;
 * @param res contains the root on exit
 * @param n_max maximum number of iterations
 * @param func a function pointer
 * @return CML_SUCCESS or CML_FAILURE if n_max reached
 */
int
cml_root_bisection(cml_function_t *func, double xmin, double xmax,
                   double epsrel, double epsabs,
                   size_t n_max, double *res);

/**
 * Find the root of a function using Newton's algorithm with the Armijo line
 * search to ensure the absolute value of the function decreases along the
 * iterations.  The descent direction at step k is given by
 *
 *     \f[ d_k = f(x_k) / f'(x_k) \f]
 *
 * Determine \f$ \alpha_k = \max\{2^{-j}, j \ge 0\} \f$ s.t.
 *
 *   \f[  f(x_k + \alpha_k d_k) \le f(x_k) (1 - \omega \alpha_k) \f]
 *
 * where in this implementation \f$ \omega = 10^{-4} \f$.
 *
 * @param cml_function_fdf_t pointer
 * @param x0 initial guess
 * @param x_eps if the relative improvement over the root is less than this value,
 * then stop;
 * @param fx_eps if |f(x)| < fx_eps * then stop;
 * @param n_max maximum number of iterations
 * @param res contains the root on exit
 * @return CML_SUCCESS or CML_FAILURE if n_max reached
 */
int
cml_root_newton(cml_function_fdf_t *func, double x0,
                double x_eps, double fx_eps,
                int n_max, double *res);

ScientificC / cmathl